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We systematically study a collection of refractive phenomena that can 
possibly occur at the interface of a two-dimensional photonic crystal, with 
the use of the wave vector diagram formalism. Cases with a single propa- 
gating beam (in the positive or the negative direction) as well BjS crises with 
birefringence were observed. We examine carefully the conditions to obtain 
a single propagating beam inside the photonic crystal lattice. Our results 
indicate, that the presence of multiple reflected beams in the medium of in- 
cidence is neither a prerequisite nor does it imply multiple refracted beams. 
We characterize our results in respect to the origin of the propagating beam 
and the nature of propagation (left-handed or not). We identified four dis- 
tinct cases that lead to a negatively refracted beam. Under these findings, 
the definition of phase velocity in a periodic medium is revisited and its phys- 
ical interpretation discussed. To determine the "rightness" of propagation, 
we propose a wedge-type experiment. We discuss the intricate details for an 
appropriate wedge design for different types of cases in triangular and square 
structures. We extend our theoretical analysis, and examine our conclusions 
as one moves from the limit of photonic crystals with high index contrast 
between the constituent dielectrics to photonic crystals with low modulation 
of the refractive index. Finally, we examine the "rightness" of propagation 
in the one-dimensional multilayer medium, and obtain conditions that are 
different from those of two-dimensional systems. 
PACS numbers: 78.20.Ci, 41.20.Jb,42.25.-p, 42.30.-d 
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I. INTRODUCTION 



Photonic crystals (PC's) are dielectric structures with two- or three- di- 
mensional periodicity. They are known as the semiconductor counterpart 
for light, since they exhibit the ability, — when engineered appropriately — 
to mold and control the propagation of EM waves. Among their unusual 
properties lies their ability to exhibit a wide variety of anomalous refractive 
effects, which recently attracted a great deal of interest, both theoretically 
[H 121 El HI E] and experimentally [UJ [7j. The observed refractive effects can 
be quite complicated, and, in most cases, the direction of the propagating 
signal cannot be interpreted with the use of a simple Snell-like formula. In 
particular, Kosaka et al. [Hj observed a large swung of angle for the refracted 
beam, for a small angle of incidence. They called this effects the "superprism 
phenomenon." 

Anomalous refractive phenomena are known in the field of optics and are 
commonly associated with anisotropy in the optical properties of the material 
(permittivity) jH]. Two propagating solutions exist, having a different disper- 
sion relation. One of them is extraordinary, i.e., non-spherical. As a result, in 
some cases two refracted beams are observed in these media, a phenomenon 
known as "birefringence" (HI E] • Numerous studies GO HH C21 EH 03] on 
diffraction gratings, essentially the one-dimensional (ID) counterpart for the 
PC structures, led to the observation of a vast variety of anomalous refracted 
effects, including "birefringence." These systems have undergone extensive 
and systematic study [TUl HD U21 EE El , based on the wave vector diagram 
formalism. This formalism was proven to be an excellent tool in explaining 
the unusual refractive properties for the ID diffraction grating system. The 
reader can find a didactic description of these diagrams and their use in Refs. 
GDI and GH- 

Despite the recent interest focused on the superrefractive effects in two- 
dimensional crystals, a systematic study is certainly lacking in the literature 
for these systems and only a few effects were studied and discussed Gl 121 El El 
El US]. So far, all studies on the two-dimensional (2D) photonic crystals are 
restricted to effects that can be explained with the study of the corresponding 
propagating modes within the first Brillouin zone (BZ). However, as we will 
demonstrate in this paper, a class of unusual propagation phenomena in PC's 
can be explained only by a careful study of all allowed propagation modes 
in all zones. Our analysis indicates that, contrary to one's intuition, the 
multiplicity of reflected beams in the incoming medium does not necessarily 



2 



imply the presence of multiple beams in the PC medium and vice versa. In 
this context, we also investigate carefully the conditions necessary to obtain 
single beam propagation inside a 2D square or triangular PC lattice. 

The anomalous refractive effects observed in both the ID grating and PC 
literature include cases where the light bends "the wrong way," i.e., is re- 
fracted negatively at the air-PC interface. Such a phenomenon was observed 
and widely discussed in the left-handed materials literature (TBI [T71 [TBI EH] ■ 
In the left-handed medium (LHM), homogeneous jT"""| or composite [T7|, the 
electric field vector E, the magnetic field vector H, and the wave vector k 
form a left-handed set of vectors. The sign of the product S ■ k — S being 
the Poynting vector — reflects the sign of the "rightness" for the system 
and is negative for the left-handed medium (LHM). It is also customary 
to refer to a left-handed propagating wave as a backwards wave [20] ■ The 
refractive index for such a medium was calculated with the use of the scat- 
tering data and was found to be unambiguously negative [JJ"]. However, the 
characterization of the left-handed or right-handed nature of propagation is a 
point somewhat overlooked both in the ID gratings and PC literature. Only 
in Ref. [5] the sign of the product S • k j2JJ, with k in the first BZ, was 
determined for a two-dimensional hexagonal PC with a finite difference time 
domain simulation (FDTD) [22] • The simulation experiment, performed on 
a wedged PC structure, is in accordance with the UCSD experiment [21"", 12"""] . 
In the latter, the negative index was experimentally verified for the tradi- 
tional composite LHM. Left-handed behavior in photonic crystals relates to 
the origin and nature of a certain propagating beam. We intend to study 
this for different cases with the wave vector diagram formalism. In any case, 
the assignment of a proper refractive index, should carry the information 
regarding the rightness of the PC medium in its sign and be consistent with 
the left-handed literature. 

Moreover, the phase velocity for an EM wave propagating in a periodic 
structure is a subject of some controversy in the literature. Yariv defined 
the phase velocity for a propagating EM wave in the ID layered medium as 
the phase velocity that corresponds to the dominant plane wave component 
|25j . Recently the phase velocity has been associated with the Bloch's crystal 
momentum k, where k is in the first BZ ["fl [""J "2"""| • Specifically, in Ref. [JJ the 
phase velocity and appropriate phase index were discussed for both limits of 
index contrast between constituent dielectrics (high and low). Considering 
this controversy, the subject of phase velocity in a periodic medium should 
be revisited. It is certainly worthwhile to reexamine the physical meaning of 
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each definition in both limits of refractive index modulation (high and low). 

In an attempt to make the study of the PC system simpler, in some cases 
the PC system was homogenized appropriately with the use of an effective 
medium theory [27] . However these theories mainly apply to the long wave- 
length limit. However, as we will show in our subsequent analysis, for some 
cases that lie in the higher bands, it is still possible to characterize the refrac- 
tive and propagation properties with an effective index n(u) under certain 
conditions. It is important to examine carefully such conditions, since the 
study of the PC can be greatly simplified. We will see that in these cases, 
both phase and energy velocities can be derived by simple formulas. 

In this work, we attempt a systematic study for the anomalous refractive 
phenomena occurring in two-dimensional PC systems. We focus on various 
cases that have substantially different origins and nature. The characteristics 
of each case are analyzed. For this purpose, we discuss the properties of 
phase and energy velocities, as well as types of propagation (left- or right- 
handed). In particular, in Sec. II we present four distinct cases of anomalous 
refractive effects, where a negatively refracted beam is present. We explain 
and analyze the origin of the refracted beam with the wave vector diagrams in 
Sec. III. Using the same formalism, we characterize birefringent phenomena 
in photonic crystals that we observed. We also discuss how these relate to 
the presence of Bragg reflections in the medium of incidence. We discuss all 
relevant properties for an EM wave propagating in the crystal. In particular, 
we define appropriately a phase velocity, calculate it numerically, and discuss 
the meaning of the associated effective phase index in Sees. IV , V and VI, 
respectively. Moreover, in Sec. VI we discuss the conditions necessary to 
obtain single beam propagation. We derive expressions for both group and 
energy velocities, and show their equality in Sec. VII. In Sec. VIII, we 
discuss the group refractive index associated with the group velocity. In Sec. 
IX we focus our discussion on the left- or right-handed nature of propagation. 
In this section, we discuss the details of the appropriate wedge experiment 
design, that unveils the sign of "rightness" for propagation inside the PC 
for different cases of triangular and square lattices. In Sec. X, we discuss 
the validity of our theoretical analysis as one moves from the limit of high- 
index modulated crystals to the limit of photonic crystals with low index 
modulation. Finally, we make a comparison between the two-dimensional PC 
medium and the ID layered medium in Sec. XI. We present our conclusions 
in Sec. XII. 
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II. ANOMALOUS REFRACTIVE PHENOMENA AT THE AIR 
PC INTERFACE 



We present in Fig. 1 four characteristic cases where a negatively refracted 
beam appears when light is incident at a PC slab interface. The cases shown 
in Fig. 1 basically outline the different possible reasons for which a negatively 
refracted beam can appear inside the photonic crystal. In the case of Fig. 
1(b), two distinct beams propagating in opposite directions (positive and 
negative) are present (birefringence). To study the various super- refractive 
effects, we employed the Finite Difference Time Domain technique (FDTD) 
[22 I2H] with Perfect Matched Layer (PML) 29J boundary conditions. We 
study various triangular PC structures of dielectric cylindrical pillars in air 
for the H (TE)-polarization case (magnetic field aligned along the cylinder's 
axis). Whenever possible, we used the value of 12.96 for the dielectric con- 
stant and r=0.35 for the radius of the rods for consistency and comparison 
with the results in Ref. [5] and the results of Notomi pQ. However sometimes 
for the purpose of isolating and observing clearly specific effects, it becomes 
necessary to employ PC structures with different parameters. The presence 
of a negatively refracted beam is clear in all four cases as seen in Figs. 1(a)- 
(d). Before we expand our analysis, we discuss the wave vector diagrams. 
Careful use of such diagrams in the PC system can always explain/determine 
the direction(s) of the refracted beam(s). Then, we will be able to comment 
on the nature/origin of each different superrefractive effect shown in Fig. 1. 



III. WAVE VECTOR DIAGRAMS AND INTERPRETATION 
OF THE FDTD RESULTS 

The wave vector diagram contains the equifrequency surfaces (EFS) for 
the photonic crystal that apply for the frequency of operation. Actually, for 
our two-dimensional system the surfaces reduce to contours. These contours 
consist of all allowed propagation modes in wave vector space that exist in 
the PC system for a certain frequency. One or multiple contours can be 
relevant for a certain frequency, depending on the number of bands corre- 
sponding to the frequency of interest. To isolate the different effects, we 
focus our study on cases with only one band corresponding to the frequency 
of interest. Thus, we have only a single EFS within the first BZ in k-space, 
closed or broken (with six-fold symmetry for the hexagonal lattice). The cor- 
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responding analogue of the EFS for the electronic case would be the Fermi 
surface. However, unlike electronic states that exist inside the crystal, an EM 
wave propagating inside the photonic crystal is excited with an EM wave in- 
cident from the outside. This implies the following: 1) the propagating state 
inside the crystal has the same frequency as the frequency of the source. 
2) the causal direction of propagation inside the crystal points away from 
the source. 3) the wave vector of the propagating wave inside the crystal is 
subject to additional restrictions, imposed by the boundary. Specifically, the 
parallel component of the wave vector is given by the following formula |3U] • 



'sir, cut 
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k\\,m = - Sin d inc + r , (1) 

C Vstr,cut 

where we consider EM waves incident from air with angle 9i nc measured from 
the surface normal, m is an integer equal to 0, ±1, ±2 and b str ,cut represents 
the periodicity of the interface. For different interface cuts and lattice ar- 
rangements, we have the following cases: 

for triangular cut along TK 
for triangular cut along TM , (Al) 
for square cut along TM 
for square cut along TX 

with a being the lattice constant. 

Basically, formula (1) is the generalization of the phase matching con- 
dition at a periodic boundary [3Tj- So, the kn conservation condition as 
expressed with Eq. (1) is an integral part of the wave vector diagram. Dia- 
grammatically, Eq. (1) can be represented by m parallel lines, all perpendic- 
ular to the interface and separated by 27i/b str , C ut- We refer to these lines as 
construction lines in accordance with existing nomenclature in the literature 
(see for example Ref. [TU]). 

The EM wave that propagates in the photonic crystal is a superposition 
of many plane waves, called Floquet-Bloch wave (FB wave) ^Hl E2] . It is 
characterized by the wave vector in the first Brillouin zone, referred to in 
the following as the fundamental wave vector We will discuss the char- 
acteristics of the FB wave in detail in Sec. IV. The intersections between 
EFS and construction lines represent the possible wave vector values for the 
FB propagating beam(s). However, not all of these intersections correspond 
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to a signal pointing away from the source (causal signal). So, for each wave 
vector intersection we need to determine the direction of the corresponding 
signal, which is nothing but the direction of the corresponding energy ve- 
locity We will show in Sec. VII that the energy velocity v e is equal to 
the group velocity v g for the photonic crystal. Now, v g is VkV- Thus, the 
geometric properties of the gradient require the propagating signal to have 
a direction normal to the EFS at a certain wave vector point in k-space and 
to point towards increasing frequencies to. Accordingly, before we proceed 
with a wave vector diagram analysis, additional knowledge for the sign of the 
product Vk^ • k is necessary. Essentially, we need to know the curvature of 
the relevant band. 

As already mentioned in the introduction, wave vector diagrams have 
been used before in the PC literature [1112113130113, but were drawn in the 
first Brillouin zone only. Although in many cases this methodology suffices, 
for a complete treatment such diagrams should include the EFS shown in 
the repeated zone scheme. A complete treatment yields all possible waves 
that can couple inside the photonic crystal. All intersections between the 
construction lines and the EFS in the repeated zone scheme should be ac- 
counted. Apparently, in our repeated zone treatment, some intersections 
may have wave vectors falling outside the first BZ. These, should be folded 
back to the first zone, in order to obtain the fundamental wave vector of the 
FB wave. The folding process involves adding or subtracting an appropriate 
reciprocal lattice vector. Actually, in our study cases, where the interface is 
chosen along a symmetry direction, only the primary construction line (the 
one with m=0) is sufficient. We refer to it, as simply construction line or 
k|| conservation line. To convince the reader, we show an example (see Fig. 
2) with all the construction lines present. However, for simplicity, in all the 
other diagrams the higher order construction lines are dropped. Summariz- 
ing, in any case we have fixed frequency (EFS contour), fixed sign Vk^ • k, 
and fixed parallel component of wave vector (construction line) . Our method- 
ology stems from the properties of the FB wave, which we discuss in more 
detail in Sec. IV. Following the folding process many points may fall onto the 
same point in the first zone. We emphasize that in this case, all points give 
rise to one beam only. We refer to these points as "equivalent" points. Only 
different wave vector points in the first zone (after the folding process) yield 
different beams, provided the corresponding signal is causal. We distinguish 
between effects that stem from the first and the higher Brillouin zones. A 
beam that originates from k points in the first BZ is the "transmitted" beam, 



7 



while we refer to beams coming from k points in the higher order zones, as 
higher order beams. 

The PC medium can also give rise to multiple reflected beams appearing 
in the medium of incidence. We can choose to determine these graphically. 
In this case, one would need to keep all construction lines determined by Eq. 
(1). Alternatively, the angle of an order m reflected beam, when the medium 
of incidence is air, is given by the following simple formula |3U| : 



e m = sin" 1 % (2) 



provided 



kim < u 2 /c 2 . (3) 

(with fc[|, m given by Eq. (1)). If condition (3) is not satisfied Vm ^ 0, then 
no additional higher order reflected beams exist. To avoid confusion with 
the higher order waves in the PC medium, we refer to the latter as higher 
order Bragg reflected beams. Notice that formula (2) also provides the angle 
for the order m outgoing beam, in the case that air succeeds the PC slab 
material. 

In Fig. 2 we show the wave vector diagram for the case of Fig. 1(a) 
drawn in the repeated zone scheme. Note that all the EFS are calculated 
with the use of the plane wave expansion method. We note that whenever we 
refer to the PWE method [3U 023 EE] for the H-polarization case, we applied 
Ho's method instead of the inverse expansion method, since the former is 
proven to show faster convergence (S3 EE] • The bold green dot-dashed line 
in Fig. 2 represents the construction line. It intersects points A, B of the 
EFS in the first zone (black circle) and points A2, B2, A3, B3 of the EFSs 
in the higher order zones. We fold points A2, B2, A3, B3 back to the first 
zone by adding G n =n G oy (see figure) (where n=-2 for points A2, B2 and 
n=+2 for A3 and B3). Notice they all fall back onto points A and B. The 
case of Fig. 2 corresponds to a band with negative curvature. Therefore, A 
has v e pointing away from the source, while B has v e pointing towards the 
source. This means that only point A contributes to a propagating beam, 
the "transmitted" beam. We indicate the propagating beam with the bold 
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orange arrow in Fig. 2, and the bold solid arrow in the FDTD simulation in 
Fig. 1(a). As we mentioned earlier, only the zeroth order construction line 
is really necessary in the determination of the propagating beams. However, 
in Fig. 2 for didactic reasons, we plotted the higher construction lines as 
well (dotted cyan vertical lines). Each is displaced in respect to the zeroth 
order one (bold green vertical line) by Gi nt = A G ox (where A = ±1,±2 
e.t.c). We plotted in the figure only the higher order construction lines with 
|A| = 1 . We account for all intersections of these lines. Subsequently, we fold 
them back to the first zone by adding/subtracting an appropriate reciprocal 
lattice vector (G=n G oy +A G ox with rail, Ail). All fall back onto 
either A or B. Evidently, only the zeroth order construction line is sufficient. 
We checked and determined this generally holds for any two-dimensional PC 
slab, provided that the interface is cut along a symmetry direction. From 
now on we draw only the zeroth order construction line. In the insert of Fig. 
2 we show a wave vector diagram for the same case, but drawn in the first 
zone only. Apparently, in this case an analysis in the first BZ gives identical 
results with an analysis with the modes drawn in the repeated zone scheme. 

In Fig. 1(b) we see that two beams coexist (birefringence). This case 
corresponds to a band with positive curvature (v e • k > 0). A wave vector 
diagram in the repeated zone scheme is seen in Fig. 3. The green bold solid 
line in the diagram is the zeroth order construction line for this case. We 
mark the intersections of this line with the modes in all the zones and fold 
back those falling outside of the first zone (by adding G= -G ox in G oy with 
ra = il). Notice, points A2, B2, A3, B3 when folded back in the first zone, 
fall onto points A', B' that are different from A and B. From the set of wave 
vector points resulting in the first zone (A, B, A', B') only B and B' give 
a signal that propagates away from the source. The respective signals are 
indicated with the bold (point B) and dotted (point B') orange vector in 
Fig. 3. They correspond to the bold and dotted black arrows in the FDTD 
simulation of Fig. 1(b). Actually, the first beam (solid arrow in Fig. 1(b)) 
can be explained with an analysis in the first zone (see insert of Fig. 3) 
and is the transmitted beam. Nevertheless, the beam that corresponds to 
the dotted arrow clearly stems from a higher order zone, i.e., a higher order 
effect. The latter would be impossible to predict with an analysis within the 
first Brillouin zone. 

In Fig. 4 we see the EFS plotted in the repeated zone scheme for the case 
of Fig. 1(c). The construction line (bold green solid line) intersects point A 
and B in the first Brillouin zone and points Al, Bl, A2, B2 in the higher 
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order zones. These fall onto points A' and B' when they are folded back to 
the first zone (with a similar process as in Fig. 3). Taking into account the 
sign of v e • k, in this case positive, from all intersections A, B, A', B' only B 
and B' yield a propagating beam (FB wave) (bold and dotted orange arrow 
in Fig. 4). The one that corresponds to point B can be determined with an 
analysis within the first BZ (see insert) and is the transmitted beam. The 
second one is a higher order effect. The origin for this higher order wave 
is similar to that of Fig. 3. However, because the EFS are anisotropic in 
this case, both beams (transmitted and higher order) are negatively refracted 
beams. These types of higher order waves, when EFS are anisotropic and 
broken, are unique to the triangular structure. They stem from the six-fold 
symmetry of the modes in the wave vector space. For square structures with 
equifrequency contours broken (4- fold symmetrical in this case) , birefringent 
effects of this kind cannot be observed. 

In Fig. 5 we show the wave vector diagrams in the repeated zone scheme 
that corresponds to the case of Fig. 1(d). The ky conservation line intersects 
several points (A2, B2, A3, B3). However all of the points, when folded back 
in the first zone, fall onto either point A' or B' (are equivalent to A' and B' 
respectively). Since v e • k >0, only the wave vector at B gives a propagating 
FB wave, with energy velocity indicated by an orange solid arrow in the 
figure. This is the sole propagating beam shown with a black dotted arrow 
in the corresponding FDTD simulation in Fig. 1(d) and is essentially a 
higher order beam. Note, an analysis within the first BZ (see insert of Fig. 
5) predicts no propagating beam in this case. A wave vector analysis in the 
repeated zone is needed. 

Our preceding analysis shows that in any general 2D PC system, the di- 
rection^) of the propagation beam(s) can always be determined with careful 
use of the wave vector diagrams in the repeated zone scheme. Notice the 
excellent agreement between the theoretical prediction — derived from the 
wave vector diagrams — and the actual FDTD simulations seen in Fig. 1. 
We do not discuss the aspect for the "rightness" of propagation in this sec- 
tion. We will see in Sec. IX that cases with v e • k < represent a backwards 
(left-handed) wave. We note at this point that only the negatively refracted 
beam in Fig. 1(a) is a backwards wave. Actually, the mechanisms that lead 
to the formation of a negatively refracted beam in the cases of Figs. 1(a) 
through (d) are distinctly different. In particular, in Fig. 1(a) a negatively 
refracted beam is formed, because the perpendicular component of the wave 
vector reverses sign when meeting the air-PC interface. In fact, this is the 
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same mechanism leading to a formation of a negatively refracted beam in the 
homogeneous negative index medium (NIM). However, this is not the case 
for Figs. 1 (b) through (d). Notice that in all these cases the perpendicular 
component of the wave vector does not reverse sign at the PC interface. For 
the case of Fig. 1(c), the negatively refracted beam , indicated with the 
solid arrow, is due to the anisotropy of the EFS in k-space. The negatively 
refracted beams in Figs. 1(b) and 1(d) stem from modes in the higher order 
zones. In other words, both are a higher order effect. The same holds for the 
second negatively refracted beam in Fig. 1(c) (dotted line). Nevertheless, 
still the nature of the higher order waves in Figs. 1(b) and 1(d) (dotted lines) 
is different. In the case of Fig. 1(b), the component of the incident wave 
vector parallel to the interface kii falls within the limits of the surface ID 
BZ. This zone extends between — TT/b str ,cut, n/b s t r ,cut, with b str ,cut being the 
interface periodicity given by expression (Al). On the contrary, in the case 
of Fig. 1(d), k|| falls outside the ID surface BZ limits. In order to distinguish 
between the two waves, we refer to them as type I and type II, respectively. 
Type II higher order waves are present because of the periodicity across the 
interface. Type I higher order waves are present because of the periodicity 
of the whole bulk photonic crystal. Hybrid higher order effects of type I and 
II can also be observed in some cases. Type II waves are generic and can be 
observed for any case, provided that frequency and/or angle of incidence is 
high enough. Type I higher order waves are particular to the specific lattice 
type, symmetry direction of interface, and EFS features. 

We emphasize that higher order waves of type I always coexist with a 
beam deriving from the first BZ (transmitted beam). Thus, whenever a type 
I wave is present, birefringence is observed. This is the case of Figs. 1(b) and 
1(c). Birefringent effects in photonic crystals were observed experimentally 
by Kosaka et al. in their study, for the frequency where the double 

branching of the beam was observed, two-band solutions exist. In other 
words, two different equifrequency contours exist within the first Brillouin 
zone for the relevant frequency. When multiple bands exist for a certain fre- 
quency, the procedure we just described in detail must be repeated for each 
separate band. Many more beams can propagate in these cases. We note 
that Born and Wolf 8 J (as well as Yariv and Yeh jSTj) adopted an effective 
medium approach to describe the ID layered medium. They found that it 
effectively behaves as a homogeneous medium with optical anisotropy. There- 
fore, the ID layered medium is capable of showing birefringent effects. They 
termed these effects as "form" birefringence to stress the fact that these orig- 



11 



inate from anisotropy on a much larger scale than the molecule. In optically 
anisotropic materials, the tensor property of the permittivity introduces two 
solutions for the dispersion relation (ordinary and extraordinary). The two 
different dispersion relations give two different equifrequency surfaces in the 
wave vector space and lead to the familiar birefringent phenomena in these 
media. In a way, we can say that multi-fringent phenomena in PC's arising 
from multiple bands appear for similar reasons as the ones in the optically 
anisotropic materials. In essence, multiple bands imply multiple dispersion 
relations for a certain frequency region, and, therefore, multiple EFS within 
the first BZ. However, all may be extraordinary (non-circular in 2D). In this 
paper, we observe in Figs. 1(b) and 1(c) birefringent effects that have a to- 
tally different origin. In these cases, there is a single band for the operation 
frequency, i.e., a single branch for the dispersion relation, and, therefore, a 
single EFS exists in the first BZ. The two beams arise because the modes, 
represented by the EFS, repeat themselves periodically in the wave vector 
k-space. This means that the periodicity of the PC lattice comes into play in 
two different ways as far as birefringent effects are concerned. It introduces 
the possibility of having multiple dispersion relations within the first BZ for 
a certain frequency region, as in the case of Kosaka et al. ^3] . On the other 
hand, the periodicity causes the modes to repeat themselves in reciprocal 
space. This is responsible for the beam doubling effects we observed in Figs. 
1(b) and 1(c). 

The existence of multiple beams in the medium of incidence and the 
medium succeeding the PC slab, as well as how these relate to the trans- 
mission properties, were previously studied [SUJ. Despite these studies, the 
relation between the existence of multiple beams in the medium of incidence 
and inside the PC has not yet been carefully examined. Incidentally, Luo et 
al. |3] state that the condition necessary to obtain single beam propagation 
inside the PC is u ^ 0.5 x 2ixc/a s , where a s is the interface period. The 
quoted condition can be rewritten as / ^ a/2b strtCUt , with b str ,cut given by 
(Al). In fact, if EM waves are incident in the PC slab from air, this con- 
dition guarantees the absence of any higher order Bragg reflected beams for 
any angle of incidence (see Appendix 1(a)). For a triangular lattice cut along 
TM, this condition becomes / ^ 0.289. However, FDTD simulation results 
that we present in the following show this condition does not guarantee 
single beam propagation inside the PC medium. We consider the case of 
dielectric rods with permittivity e = 60. and radius r=0.37 a. This is a case 
qualitatively similar to that of Fig. 1(b) (Fig. 3), but with a much lower 
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relevant frequency (/ = 0.275). We stress that one band only corresponds 
to this operation frequency. Evidently / = 0.275 is below the quoted limit, 
which means no higher order Bragg reflected beams exist for any angle of 
incidence. Indeed, in Fig. 6 we observed only one reflected beam. Notice, 
however, the clear presence of two propagating beams (solid and dotted ar- 
row). The second beam (dotted arrow) is a higher order wave of type I, like 
the one shown in Fig. 1(b). 

A single beam propagation condition cannot be derived always in a sim- 
ple manner and one should, in general, carefully examine the wave vector 
diagrams in the repeated zone scheme. When after the folding process, only 
a sole wave vector in the first BZ gives a causal FB wave, only then do we 
have single beam propagation. Consequently, the presence of only a single 
reflected beam is neither a prerequisite nor does it guarantee the presence 
of a single beam coupling into the PC medium. Also, note in Fig. 1(d) 
the clear presence of a higher order reflected beam, while there is only one 
propagating beam. For certain simple cases of isotropic EFS, we will discuss 
the conditions for single beam propagation in Sec. VI. 



IV. FLOQUET BLOCH WAVE AND PHASE VELOCITY 

Consider the magnetic field of an H-polarized wave inside a two-dimensional 
periodic photonic crystal structure for the case of H (TE) -polarization. 

H(r, t) = ~^=^ r V H G (k, ^ k )e' G - r e-^. k * z. (4) 
V A ws q 

A\ys is the area of the Wigner-Seitz cell, z is the unit vector in the direction 
out of the plane of periodicity (i.e., the direction of the cylindrical rods) and 
G is a reciprocal lattice vector. The coefficients are determined from the 
eigenvalue equations obtained from the PWE method jnHESEHl- Apparently 
the above expression for the field satisfies the Floquet-Bloch theorem 
for a periodic medium. A field that propagates according to expression (4) 
is known as a Floquet-Bloch wave |3*2"l ITUj with k lying in the first zone. 
Any attempt to express the propagation solution in terms of wave vectors 
k lying outside the first BZ results in an expression equivalent to Eq. (4) 
(see Appendix 1(b)). So, the wave vector chosen in the first zone is what 
characterizes a propagating FB wave. We call this the fundamental wave 



13 



vector [TOJ ESI- Hence the term, "equivalent points," describing points in 
k-space separated by a reciprocal lattice vector. This property of the FB 
wave explains the general recipe that we followed in the preceding section 
to determine the propagating waves. Clearly in the 2D periodic system, all 
plane wave components contributing to the FB wave with fundamental wave 
vector k (expression (4)) propagate together, not separately, with a common 
energy velocity, v e . Note, no individual plane wave components serve as a 
separate solution of Maxwell's equations. As a result we do not see clear 
phase fronts, but rather have phase-like fronts with a "wiggly" profile (see 
for example the FDTD simulation presented in Fig. 1). This type of profile 
manifests a strong plane wave component mixing PQ, present in PC crystals 
with high refractive index modulation. 

The questions of how one should approach the subject of defining a phase 
velocity for the FB wave is still unanswered. What would really be the 
physical meaning for such definition. Yariv [25] defined a phase velocity 
for a ID periodic system (see Appendix 1(c)). Equivalently, for the two 
dimensional system this definition would be 

o 

with K Q = k + G being the plane wave component that has the larger am- 
plitude Hc o in expression (4). In other words, it is the wave vector of the 
predominant plane wave component. Many PC studies [211 EH EE] focused 
on the long wave length limit, where such a definition would be appropriate 
|25j . Nonetheless, the interesting refractive behavior of the photonic crystal 
reveals itself in the higher bands. Unavoidably, the subject of phase velocity 
in the photonic crystal requires some rethinking. As a matter of fact, Notomi 
[T] , as well as Kosaka [2H] , defined a phase index that corresponds to the fun- 
damental wave vector within the first BZ, of the FB wave. Accordingly, the 
phase velocity would be: 

v P = ^k, (6) 

where k is in the first BZ. 

There is an apparent contradiction between these two definitions as given 
by expressions (5) and (6), respectively. To investigate for the physical mean- 
ing of the phase velocity in a periodic medium, we will study numerically the 
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field patterns of the propagating wave in the next section. 



V. NUMERICAL DETERMINATION OF THE PHASE VE- 
LOCITY 

We consider two cases of almost isotropic EFS. In both cases the fun- 
damental wave vector of the propagating FB wave is chosen to lie along a 
symmetry direction. We choose again the structure of Notomi (i.e., the same 
structure as in Fig. 1(a)) and two frequencies a) cji=0.58 2nc/a lying in a 
band with negative curvature and b) co>2=0.48 2irc/a lying in a band with 
positive curvature. In order to extract information about the phase veloc- 
ity in the system, we need to compare the time-independent fields at various 
points along the propagation direction. This methodology is analogous to the 
one followed by Ziolkowski and Heyman j^H], where the negative phase index 
was numerically confirmed for a homogeneous slab material with e = — 1 and 

For this purpose, we consider a pulsed signal f(t) cos(u;oi) with 



/(*) 



if t < tl 



2 +{t-n) 2 
(t-fi) 2 



a 2 + (t-t3) 

if t > f 3 



if tl < t < t2 . (A2) 
if t2 < t < t3 



The parameters are chosen to give a broad pulsed signal in time, with 
a small 5u around the operation frequency u> . For operation frequency u 1 
the parameters are a=21.9 T, tl = 19.6 T, t2=499.7 T, and t3=979.9 T. For 
operation frequency co>2 they are a=18.1 T, tl=16.2 T, t2=413.6 T, and 
t3=810.9 T. The period T of the EM wave is different for the two cases and 
the parameters are chosen to correspond to the same actual time. The pulse 
is launched normally to the photonic crystal structure, with an interface 
cut along the TK symmetry direction. We monitor the magnetic field H 
for each time step for a long time, for certain points along the propagation 
direction. We refer to these points as detector or sampling points. The 
Fourier transform of the time series H(y,t), where y the coordinate of a 
detector point, yields the corresponding amplitude H(y, uj). 

Before we proceed with our analysis, we should mention that in order to 
make any assessment regarding the phase velocity the field patterns inside 
the slab should be as close as possible to the infinite system patterns, given 
by the FB wave expressed in (4). For this purpose, our structure should 
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emulate a semi-infinite PC slab. Any wave that couples into the slab will 
undergo multi-reflections between the two interfaces. In order to achieve our 
goal, i.e., a structure acting similar to a semi- infinite slab, we must somehow 
eliminate or reduce substantially the amplitude of any reflected beam orig- 
inating from the second interface. For this purpose, we consider a periodic 
structure consisting of 30 sites along the lateral direction and 100 rows along 
the propagation direction. Of these 100 rows, 30 rows consist of rods with 
dielectric constant 12.96 embedded in air. This part of the structure is the 
area of concentration. We take detector points that monitor the field as a 
function of time in this area. In the remaining 70 layers, we introduce absorp- 
tion (both in the sites and in the background) so that the field is attenuated 
before exiting the crystal. We introduce electric and magnetic conductivity 
(<T e and a m , respectively) to each numerical cell, so that the impedance of the 
each grid cell in the absorptive layer would be the same as the corresponding 
one in the non-absorptive layer. It is equal to ^J-^f—, where eo and fiQ are 

the vacuum permittivity and permeability, and q j is the relative permittiv- 
ity of the 2D grid cell located at the point with grid coordinates (i, j). This 
is possible when the conductivities that we introduced follow the relations: 
a e = e it j(T and a m = ^-a , where a is a conductivity parameter. If the 
parameter (Jq is chosen very low (lO~ 5 Q~ 1 m~ 1 ), the reflections of the beam 
when entering the absorptive layer are low. In order to make a rough estimate 
of the EM wave energy that gets reflected back to the non-absorptive layer 
where we monitor the fields, we look at the attenuation profile of the fields 
inside the absorptive layer. In addition, reflections can occur when the EM 
wave meets the absorptive boundary. To check these, we considered oblique 
incidence. Overall, we find the amplitude of the field that gets reflected back 
into the non-absorptive layer is about 10% of the amptitude of the refracted 
EM wave. 

The set of points that serve as detectors are chosen normal to the surface, 
which coincides with the propagation direction y, chosen along the TM sym- 
metry direction. Their respective locations are a distance b = \/3a apart, 
essentially the periodicity for the propagation direction. We take the Fourier 
transform of the time series representing the evolution of the magnetic field at 
a certain point jUJj. Afterwards, we calculate the ratio's H(u, di+i)/ H(u>, di). 
uj represents the frequency of the input pulse train and di the location of the 
i th detector point. Since the distance between the detectors is one period 
along the propagation direction, this ratio should be equal to exp(ikb), with 
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k restricted in the first BZ. Therefore, by studying the field patterns, we 
can extract information about the phase velocity, defined in accordance with 
Notomi pQ (expression (6)). 

We calculate the field ratio at adjacent detector points and extract the 
wave vector from the following formula 



h = -Jn^^, (7) 
ib H(u>, di) 

where i stands for the i th detector point. Notice that the Fourier transformed 
fields are complex and therefore the logarithmic function is complex and 
multivalued. We record all possible values with the real parts falling inside 
the first BZ. Taking the average for the various detector points for the case of 
to = ul, we find that two possible solutions exist for k: 1) k=(1.47+0.005i) ± 
(0.01 + 0.008i) a- 1 and 2) k=(-2.16 + 0.005z) ± (0.01 + 0.008z) a' 1 . In order 
to choose the correct solution, we further study the field patterns. We also 
consider the ratio of two observation points located around the middle of 
the 30-cell PC layer, separated by Ay = 6/31. The field ratio determined 
by the FDTD simulation is 3.05 — 0.63i. Now we calculate the same ratio 
theoretically with the PWE expansion method for the two possible 

wave vector solutions determined from Eq. (7). Using solution (1) for the 
wave vector (real part only) in PWE we obtain a field ratio of 4.3 + 2.3i. 
For the second solution, we obtain a field ratio of 2.98 — 0.521 Apparently 
only solution (2) gives a ratio that agrees well with the FDTD results. We 
note that in order to eliminate any discrepancies resulting merely from the 
discretization, we used in the PWE the actual numerical dielectric grid used 
in the FDTD. However, there is still a small discrepancy between the PWE 
field ratio and numerical FDTD ratio that may stem from a combination of 
the following — angle span of incident source, Fourier transform zero padding 
errors, reflections from absorptive boundary layer, etc. This is the same 
reason for which a small imaginary part is present in the wave vector value. 
Also note there exists some ambiguity associated with the exact location of 
the first interface, since we are dealing with a periodic medium |31j . In the 
boundary layer the field values may deviate from the values given by the FB 
wave expression (Eq. (4)). Therefore, we place the first detector point at the 
center of the second row of cylinders and assign d\ — 0. 

In the analysis above we used the FDTD field patterns in the "semi- 
infinite" slab to determine that the wave vector inside the photonic crystal 
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for the case with frequency u = 0.58 x 2irc/a is k = —2.16 a -1 . This value is 
in good agreement with the corresponding value from an EFS analysis (k = 
—2.48 a -1 ) [41. Expression (6), with the use of the wave vector obtained 
from the field pattern analysis, gives v p ~ — 1.69 c y (c being the velocity 
of light). Following the same procedure but for the case with frequency 
u>2 = 0.48 x 27Tc/a, we calculate a phase velocity v p ~ 2.23 cy. In both 
cases, we determined the magnitude and the sign of the phase velocity. The 
field pattern analysis for a semi-infinite slab gives a negative phase velocity 
(i.e., opposite to the propagation direction) for the case with frequency = 
0.58 x 2nc/a). This result implies that the "rightness" of the propagating 
beam is negative in this case. So, a field pattern analysis with the FDTD 
method for a semi-infinite slab confirms the results we obtained from the 
wedge simulation experiment [S]. We will discuss more about the "rightness" 
of propagation in Sec. IX. 

In order to visualize the physical meaning for the phase velocity defined 
in Eq. (6) we plot the imaginary part of the magnetic field H{u) for various 
detector points. We show the results for both cases with frequencies oji and u 2 
in Figs. 7(a) and 7(b), respectively (open circles). The solid line is A sin ky 
where k is the wave vector as determined from the field pattern analysis 
above and y is the distance from the first detector point. Also, we choose 
an additional set of detector points, closely spaced and around the middle of 
the 30-row PC layer. We indicate the imaginary part of the corresponding 
Fourier transformed field of these points with stars in Fig. 7. We see that 
the solid sinusoidal line passes closely to the field values corresponding to the 
first set of detector points (circles). However, the field values for the second 
set of detector points (stars) deviate substantially from the sinusoidal line 
and show very high variations. 

From the FB wave expression we obtain (see Appendix 1(d)): 

< H(r = R) >= e^-*-"*) < H(r = 0) >, (8) 

where R is a Bravais lattice vector and <> the spatial average within the 
unit cell. All detectors points of the first set are Bravais lattice vectors along 
the TM symmetry direction. Our numerical results shown in Fig. 7 and 
expression (8) suggest that the phase velocity describes how fast the phase 
of the EM wave travels from period to period in the PC lattice. However, 
information of how fast the phase travels between adjacent points cannot be 
determined. Thus, it is clear why it is necessary to fold the wave vector in 
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the first zone and define the phase velocity as in expression (6). We will 
return to the same subject and the appropriateness of definition (5) when we 
discuss photonic crystals with low index modulation in Sec. X. 

VI. EFFECTIVE PHASE INDEX 

It is desirable to define an effective phase index that is correlated with 
the phase velocity as defined in Sec. IV with Eq. (6). Correspondingly, 



The sign of \n p \ is chosen to reflect the behavior left- or right-handed of 
the PC system jH] and in accordance with the left-handed literature [TB] . 
This definition for the index is consistent with the analysis by Notomi p. 
However, we have seen in Sees. II and III that the photonic crystal system 
can be quite complicated and various higher order effects may arise under 
certain conditions. One must bear in mind all these effects when interpreting 
the effective index for the photonic crystal system. Next, we analyze what 
this index represents, as well as what properties can be inferred from this 
index. 

Let's consider a Snell-like formula, 



sm9 inc = n p sm9 ref , (10) 

for EM waves incident from air, into a PC medium with phase index n p (in 
general depends on the refracted angle 9 re f). Sometimes, it is assumed that 
Snell's formula gives the direction of the propagating wave. This is not true, 
because, in general, in the photonic crystal the direction of the refracted 
wave vector and the direction of the propagating signal do not coincide. We 
discussed this in Sec. III. Note again, that the direction of propagation is 
always the direction of the energy velocity v e . Nevertheless, provided that 
certain conditions apply, there will be only a single refracted beam in the 
photonic crystal, which propagates with an angle given by Snell's formula 
(Eq. (10)). These conditions are the following: 

1) Interface cut along a symmetry direction of the crystal. 

2) Almost isotropic equifrequency contours. 

3) k inc || falling between -n/b str)CUt ,...,ii/b str:CUt 
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4) | Tip | < l/(2C7 atr , cut /) 

with b str , cut given by (Al) and 



C, 



str,cut 



1 for triangular cut along TK 

y/3 for triangular cut along TM . 

V2 for square cut along TM 

for square cut along TX 



(A3) 



Cases with slanted interfaces are more complicated to analyze, and higher 
order waves are more likely to occur, hence the first condition. The second 
condition guarantees that v e and k are about coaxial. This implies that the 
angle derived from Snell's formula represents the propagating angle. The 
third condition guarantees that the wave is not a higher order wave (specif- 
ically type II as described in Sec III). If the wave is a higher order wave of 
type II and the EFS are isotropic, it will propagate with an angle that may be 
opposite in sign to the sign of the effective index. Finally, the last condition, 
in combination with the first and second condition, guarantees the absence 
of higher order waves of type I for any angle of incidence, i.e., it guarantees 
single beam propagation. However, if we desire a single reflected beam as 
well, then the condition (see Appendix 1(a)), 



with / < a/b str ,cut, should also be observed. Notice that if / exceeds the 
value of a/b st r,cut then higher order reflected beams occur for any angle of 
incidence. 

Even in the absence of condition three, if the rest of the conditions are 
valid we still obtain single beam propagation. In this modified Snell- 

like formula can be used to determine the single refracted beam, 



Oinc < Olim — s i n 



a 



1), 



(11) 




sin 9i 



incp 



n p sin 6 re f, 



(12) 



where 




(13) 



and m chosen so that | sin6» inc - ma/(fb str>cut )\ < min(l,a/(2b strjCUt f)). 
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To summarize the purpose of defining an effective index is that it gives 
qualitative and/or quantitative insight into the PC properties such as the 
magnitude and direction of the wave vector, the "rightness" of the medium 
conveyed in the sign of the index and, under certain conditions the energy 
velocity However, caution should be taken by the use of such a phase index. 
It does not contain information about higher order beams that can couple 
inside the crystal or in the air medium (reflected beams) or both. A wave 
vector diagram type of analysis always offers a more complete treatment 
for the system. We also alert the reader that in no way should this phase 
index be used in Fresnel type formulas || to determine the transmission and 
reflection coefficients of an EM wave incident on the PC structure. 



VII. ENERGY AND GROUP VELOCITY 

We have seen in Ref. [5] that the sign of the product v e ■ k serves as a 
theoretical prediction for the sign of the "rightness" for the PC system. Such 
theoretical predictions agree well with the FDTD wedge simulation results 
jS]. We used the equality between the energy velocity and the group velocity 
to easily identify the frequency regions with negative "rightness" 0. These 
would be the regions that correspond to a band with negative curvature. 
In Ref. jSO] the equality between group and energy velocity is shown for 
3D periodic dielectric structures. It is important to show that such equality 
holds in our 2D photonic crystal as well. In order to show this, we derive 
expressions for both the group velocity v g and energy velocity v e . 

In the following we derive an expression for the energy velocity v e for the 
PC system for the H-polarization case. We start with the definition of the 
energy velocity, 

Ve = (14) 

where the brackets <> refer to the spatial average within the unit cell of 
the time averaged quantities. S is the Poynting vector and U is the energy 
density. Thus, we calculate the spatial average of the time averaged quantities 
for the Poynting vector S and the energy density U. 
The Poynting vector is defined as 

S = ^-E r x H r , (15) 
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where E r is the real electric field lying in the plane of periodicity and H r is 
the real magnetic field along the cylinder axis. The time averaged Poynting 
vector, if E and H are the respective complex quantities of the fields, is given 
then by the expression, 



with 



and 



where 



and 



S = — E x H*, (16) 
8tt K ' 



H(r,t) = e i ( k - r -^v n , k (17) 
E(r,t) = e l ( k - r -^u n , k , (18) 



v„ k ^=^H G (k,u; n , k y Gr (19) 



'A 



u n , k = -±= E G (k, oo nM ) e iG r . (20) 

Evidently the fields E and H satisfy Bloch's theorem. Note that Hq = 
Hq, z. The eigenvectors Eg and Hg are determined from the PWE method 
mSEaillH!. From (16), (17) and (18) we have 

S = ^S n , k , (21) 

with 

S„, k = u„ jk x V ; k . (22) 
Maxwell's equation for the magnetic field is 

V x H = ^e(r) E. (23) 

By substituting Eqs. (17) and (18) into Eq. (23) and by taking the cross 
product with v* k , we obtain 

(• 

S n ,k = — rr(~ iv n k x V x v n . k + < k xkx v n . k ). (24) 
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k is the fundamental wave vector of the FB wave lying in the first BZ. Then 
using expressions (19) and (24) we get 



S n ,k = — p — V (k + G)?7 G »i/ G (k,w ni k)x 



wsc „ 

G,G ,G 



H G ,(k,tu nM )e« G + G - G > r , (25) 

where r] G " are the Fourier expansion coefficients for the inverse of the dielec- 
tric function. 

Now, taking the spatial average of S n k and using the delta function ex- 
pression (see Eq.(80) in appendix II), we calculate the spatial average of the 
time averaged Poynting vector. So 



c 2 

< S >= V (k + G) 77 < vif G (k,Wn,k) H G+G ,(k,u nM ). (26) 

"' k G^G' 

For the time averaged energy density we have 

U = — (e(r)i?e(EE*) + ReCHH*)). (27) 
16n 

Taking the spatial average of the timed average quantity given in Eq. (13) 
and using the normalization conditions, 

^2 c g-g' e g • E G / = 1 (28) 

G,G' 

and 

J2H G H G = 1. (29) 

G 

We find for the spatial average of the timed averaged energy density that 

< U >= JL. (30) 
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Finally, after index manipulation, 



(k + Gi)77gi,g 2 H Gl (k,u; n;k ) H G2 (k,u n;k ). 

Gi ,G2 



(31) 



We also calculated the group velocity for this case with the use of the 
k • p (1211131111 perturbation method (see Appendix II). We found that the 
group velocity is given by the same expression. Thus, we showed the equality 
between the energy velocity and the group velocity for our 2D system. In 
fact, we also followed the same methodology for the case with the electric 
field aligned along the cylindrical rods (E- or TM-polarization). In this case 
too, we confirmed the equality between the group and energy velocity as 
well. The energy velocity for the E-polarization case is given by the following 
expression 



Evidently, the energy velocity for our 2D photonic crystal can be calculated 
with the use of the formulas (31) and (32) for the H and E-polarization cases 
respectively. One needs to know the fundamental wave vector, the index of 
the band of interest, and the FB wave coefficients (H G and E G , respectively). 
The latter are determined easily from the PWE method |3*H l3*o1 l3~o]. 



VIII. GROUP REFRACTIVE INDEX 

A group index n g can be defined jHj in accordance with traditional waveg- 
uide and optical fiber literature [IHj 



For a PC structure with the same parameters as the one in Fig. 1(a), we 
calculated the group velocity vector for the 5 th band (band with negative 
curvature) and the 4 th band (band with positive curvature) for a range of 




(32) 



c 



(33) 



(34) 
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frequencies where the EFS contours are "almost" isotropic. We used the k • p 
perturbation method |%2*1 03| E] result given in expression (82) of App. II 
a (which is equivalent to expression (31)). The results for the magnitude of 
the group velocity are shown in Fig. 8 for both bands and both symmetry 
directions (TM and TK). Notice that the closer to the band edge, the better 
the agreement between the group velocities for the two symmetry directions. 
This is expected, since the degree of anisotropy increases as one moves away 
from the band edge. Alternatively, we can consider the PC as an isotropic 
system with effective dispersive phase index n p (u). In this case, 

K\ = A, (35) 



\ n 9\ 



with, 



n g = \ n p\ + (36) 

We also show in Fig. 8 the results obtained from formulas (35)-(36) for 
comparison. Because of the small anisotropy in the EFS shape, we use for 
n p the average value of the two symmetry directions. For both bands the 
results given from Eq. (36) are shown as a solid line with circles. We see 
that the results from expression (36) are in very good agreement with the 
corresponding ones calculated from the k • p perturbation method. This is 
especially true for frequencies very close to the band edge. So, in the cases 
that conditions 1-4 of Sec. VI are satisfied, expression (35) (with the use 
of (36)) provides a good estimate for the group/energy velocity of the single 
transmitted beam. 

The sign of the group index manifests the sign of refraction at the air-PC 
interface. As in the case of the phase index n p though, caution must be 
exercised with the use and interpretations of the group index. We stress that 
the sign n g relates only to the sign of refraction for the transmitted beam. 
It does not contain any information for higher order waves of any of the two 
types discussed in Sec. III. 



IX. THE "RIGHTNESS' OF THE PC SYSTEM: DESIGNING 
THE WEDGE-TYPE EXPERIMENT 
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We observed different negative refraction effects shown in Fig. 1 in Sec. 
I. It is important to be able to characterize the nature of propagation (left- 
handed or not) for the finite PC structure. In Ref. [Sj we found that the pres- 
ence of a negatively refracted beam does not necessarily imply left-handed 
behavior. We confirmed |5] that indeed the sign of the "rightness" follows the 
theoretical prediction for the sign of n p from the band structure. A wedge 
type of experiment that can determine unambiguously the PC's "rightness" 
can be designed in most cases. However, in Sec. Ill we observed differ- 
ent higher order effects that can potentially complicate the interpretation of 
such an experiment. It is important to take into consideration the symmetry 
properties of the crystal, as they reveal themselves in the cuts of the inter- 
faces and in the eigenmodes in k-space for the frequency of operation. In 
the following, we go over the intricate details of the wedge design. We start 
our analysis with the hexagonal PC, where we identify two general classes: 
those with "almost" isotropic EFS and those with anisotropic EFS (broken 
curves with 6-fold symmetry) in k-space. Note that limiting cases between 
the two classes mentioned above exist. However, their respective frequency 
range is generally quite small and we will not treat such cases. It becomes 
evident that some a priori general knowledge for the system is necessary be- 
fore designing and/or interpreting the wedge simulation/experiment. This 
knowledge can be extracted from the band structure. One locates "almost 
isotropic cases" at the band edge, where the band structure is bell-shape like. 
Away from the band edge, anisotropic EFS are expected, which lie at the 
edge of the Brillouin zone. It is necessary to know a priori in what general 
category the shape of the EFS falls into. In fact, in some cases an estimate 
for the magnitude of the transmitted wave vector may be needed as well. 
This knowledge can be obtained from the plane wave expansion method for 
the infinite system. We stress though, that in any case the "rightness" for 
the PC from such a simulation/experiment is determined independently and 
no information regarding this quantity is borrowed from the infinite system 
analysis. 

We first consider with an almost isotropic EFS for the hexagonal 

lattices. An example is the case in Fig. 3 of Ref. [5] which corresponds 
to the system shown in Fig. 1(a) of this paper. Another example is the 
case of Fig. 1(b), that we choose to describe in this section. In Fig. 1(b) 
we observed the coexistence of a negatively and a positively refracted beam 
when the wave refracts on the PC interface. The EFS for such a case in the 
extended zone are shown again in Fig. 9. The direction for both interfaces, 
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first and wedged, are indicated with the turquoise lines. We chose TK as the 
symmetry direction for both the interfaces. From now on we will refer to a 
certain wedge design as (syml)-(sym2) where syml is the symmetry direction 
of the first interface, while sym2 is the symmetry directions of the second 
interface. So, in this case we consider a TK-TK wedge design. Notice that 
the wedge separates the space into three different areas. The area where the 
fields come from (area 1), the area inside the wedge (area 2), and the area 
after the fields experience scattering (area 3). In Fig. 9 we draw the trans- 
mitted vector inside area 2 as it would be for a system with v e ■ k positive, 
which, by the way, is the theoretical prediction. The light is normal to the 
first interface. According to the analysis in Sec. Ill, it is fairly obvious that 
there will be only one transmitted wave vector (indicated in Fig. 9 with the 
blue arrow). Naturally, in order to determine the outgoing beam we must 
apply the ky conservation at the wedged interface. We also checked in this 
case, as we did in Fig. 2, that only the zeroth order construction line (blue 
bold line in Fig. 9) suffices. We see this line intersects points PI, P2 in the 
first zone and points El, E2, E6 in the higher order zone. However, all 
points (El, E2, E6), when folded back to the first zone, fall onto either PI 
or P2. In other words, they are equivalent to points PI and P2, respectively. 
Consequently, they do not give rise to additional ky values that would cause 
additional beams^ in area 3. This is true for any TK-TK design, provided 
that | rip | < 1/(2/). We note this condition is satisfied by the vast majority 
of cases with isotropic EFS. Our subsequent analysis concerns these cases. 
The fact that all intersections El, E2, E6 are equivalent points comes 
as a virtue of the chosen symmetry direction for the wedged interface (TK). 
We bring to the readers attention that this is not, in general, the case for 
any design. In particular, TK-TM or TM-TM designs and for typical cases 
in this category, the k\\ construction line on the wedged interface intersects 
points in the higher order zone not equivalent to the intersections in the first 
zone. However, we notice these intersections, when folded back to the first 
zone, yield a k\\ value equal to k t cos9d es + 27rm/6 sirjCUt with m equal to -1 
and b str ,cut given by (Al). The angle $d es is different for different designs. It 
has the value of 60° for the TK - TM or TM - TK designs, and 6 des = 30° 
for the TK — TK and TM — TM designs. Consequently, in all forementioned 
designs, all possible ky values are predicted by the Bragg formula. Accord- 
ingly, their outgoing angle is 
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e, 



= tan 



-i 



k t cos 6 des + 




(37) 



out,m 




for m that satisfy the condition, 




(k t cos9 des + - f > 0. 



str.cut 



(38) 



We alert the reader that only for a design with interfaces along symmetry 
directions can such simple formulas be implemented. If a design is consid- 
ered with either one or both interfaces cut along a direction that does not 
coincide with a symmetry direction, then many more outgoing beams should 
be expected. 

Notice that k t is positive when parallel to the +y-direction and negative 
when anti-parallel to the +y direction. Apparently, from formula (37), the 
sign of the zeroth order outgoing beam coincides with the sign of k t . Conse- 
quently, it is the position of the zeroth order outgoing beam that determines 
the "rightness" of the PC. Since the presence of more than one beam in area 
3 may complicate the interpretation of the results, it is desirable to work 
with a design that can avoid these. In fact, for designs YK-YK and TM-TK, 
a single beam in area 3 will be present for frequencies below /o ~ 0.56. From 
these two designs it is better to choose the TK-YK. In higher bands it is very 
likely a part of the energy to get reflected back inside the PC. The result of 
these consequent multi-reflections may introduce an additional beam. If we 
choose a 60° wedge, the multi-reflected beam will always be along the normal 
to the wedge. Therefore, it is easily identifiable and should be ignored when 
present. This is the reason that the YK — YK should be used to study cases 
in this category. We performed a FDTD simulation using this design for the 
system of Fig. 9. The results are shown in Fig. 10. We observe one outgoing 
beam in the positive hemisphere. Thus, the theoretical prediction that the 
PC is "right-handed" is confirmed. The second beam around the normal is 
just the result of multi-reflections. As we mentioned above, if the frequency 
exceeds 0.56, additional beam(s) may be present in area 3. Very high fre- 
quencies for which multiple m's can satisfy Eq. (38) should be avoided in 
experimental studies. We suggest, as an upper frequency limit, the value 
of / = 0.75. For frequencies below this limit, when using the wedge design 
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indicated above and typical cases in this category, the first order Bragg wave 
has an outgoing angle still quite different in magnitude from the th order 
outgoing beam. Therefore, we can clearly identify the position of the latter 
in area 3. 

We proceed in our discussion with cases that have anisotropic EFS, i.e., 
broken curves with 6-fold symmetry in k-space. Cases that fall in this cate- 
gory are more complicated and extra caution must be exercised in designing 
and interpreting the wedge type simulation results. An example for these 
cases, is that of Fig. 4 in Ref. [5], which is the same case as the one shown in 
Fig. 1(c) of Sec. II of this manuscript. We choose to analyze another example 
belonging to this category. We consider rods with dielectric constant 12.96 
and radius 0.30a in air, magnetic field along the cylinders (H-polarization), 
and operation frequency of / = 0.50. Note, these cases have propagation 
modes only along TK. One can see in Fig. 4, for example, that a horizontal 
line across TM does not intersect any modes, neither in the first nor in the 
higher zones. This leaves us with two possible designs to work with: TM- 
TM and TM- TK design. Between the two we choose the first one, for rea- 
sons that we discuss later in this section. With the help of a wave vector 
analysis, it is easy to see that when the first interface is crossed, three differ- 
ent beams, having three different wave vectors, couple into the PC (area 2). 
One is the transmitted k t (blue bold arrow) in Fig. 11, while the other two 
are higher order beams (k dl , k d2 ). We draw the wave vectors in Figs. 11(a) 
as they would be for a "right-handed" PC and in 11(b) as they would be for 
a "left-handed" PC. At the wedged interface 12, we draw a line representing 
k|| for each of them. The blue line represents the kn value of k t , and the 
green dotted lines represent the kn values of kdi and kd2- However, in order 
to determine the outgoing angles in area 3, corresponding to the the three 
wave vectors in area 2, a simple phase matching condition represented by 
the three lines in Fig. 11, is not sufficient. To predict all possible outgoing 
beams for each of the three kii values in area 2 (kty, k^iii, kdaii), we must 
perform an analysis in the repeated zone scheme, as we did in Fig. 9. In 
fact, we should include all construction lines for each kii value to make sure 
that no possible outgoing beam is omitted. This is quite an elaborate process 
and we will not go over all the details. From this process, we observed, as a 
virtue of the symmetry the possible ky values at the interface 12 (/c||, mj i), are 
given by the Bragg formulas corresponding to each of the three wave vectors 
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(k t , k dl , k d2 ) . Accordingly, 

k\\,m,i = + 7 j (39) 

Ostr,cut 

where the index i, represents the three possible wave vectors k t kdi and kd2, 
and b str ,cut = V^a. Again, we restrict ourselves, at the lower frequencies, 
where at most the first order Bragg waves (for m=-l) can couple. An inter- 
esting observation is the first order Bragg wave corresponding to k d i has ky 
value equal to k d2 y and vice versa. In fact, this result is due to the hexago- 
nal lattice symmetry. So, in Fig. 11 we need only to add one more ky value 
(blue dotted line), that corresponds to the first order Bragg wave related to 
k t ||. The red circle in Fig. 11 represents the EFS in air. We indicate in 
Fig. 11 all possible outgoing beams for both cases of "right-handed" in (a) 
and "left-handed" in (b) photonic crystal. The bold solid and dotted blue 
arrows, and the bold dotted green arrows represent the directions of all four 
outgoing beams. Notice that the outgoing angles in each side of the normal 
to interface 12 are close. In addition, their values are quite different than the 
respective one in the opposite hemisphere. This is generally true for typical 
cases in this category, if we restrict ourselves to frequencies below ~ 0.60. 
Therefore, cases with frequencies exceeding ~ 0.60 should be avoided. Figure 
11 indicates, that the position of the larger in magnitude angles determines 
the "rightness" of the PC system. In Fig. 12 we show the corresponding sim- 
ulation for the system we chose as an example. We observe that the outgoing 
beam, with the larger in magnitude angle, lies in the negative hemisphere. 
This means that the PC is "left-handed" in this case. This agrees with the 
sign of v e ■ k obtained from the band structure. Since anisotropic modes 
reside at the edge of the Brillouin zone, the magnitude of k t will be quite 
large. Therefore, if the frequency is low, k t y may not yield an outgoing beam 
(total internal reflection). Cases with frequencies below ~ 0.50 should also 
be avoided for study. These give total internal reflection even for small com- 
paratively k t s. If one draws Fig. 11 with a 30° wedge (YM-YK design), 
we would see that we obtain outgoing beams with comparable angles in the 
two hemispheres. Thus, interpretation of the results is vague for this design, 
which led us to use the YM-YM design. 

We focused our analysis above for the hexagonal lattices. In the follow- 
ing paragraph we discuss a corresponding appropriate wedge design for the 
square structures. Our criteria are the same as in the analysis of the hexago- 
nal structures. Then, given the choice, one should always prefer the interface 
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cut for the wedged interface with the smaller periodicity, as Bragg waves (see 
formula 37) will begin to appear in such cases for larger frequencies in area 
3. Taking this into consideration and studying the modes in the extended 
zone for the square lattice, we have determined that for both classes of cases 
(isotropic and anisotropic) the appropriate wedge design is TM-TX. We 
note that the anisotropic cases for the square lattice have four-fold symme- 
try instead [3]. The four-fold symmetry does not create higher order beams 
of type I, so only one — not three as in the hexagonal case — wave vector 
couples into area 2. In addition, the wedged interface is TX in both cases 
and has periodicity a (lattice constant). For the anisotropic cases and square 
lattice PC wedge, possibly we can have only one outgoing beam in area 3. 
Typical cases may also suffer from total internal reflection, if / smaller than 
~ 0.45. We note there may be "almost" isotropic cases that suffer from in- 
ternal reflection. These will be cases that have an effective phase index with 
magnitude larger than 2 / for a hexagonal TK-YK PC wedge and larger 
than 2/v^ for a square TM-TX PC wedge. However, both of these cases 
are very rare (especially in the higher bands), since "isotropic" modes do not 
reside close to the edges of the Brillouin zone. 

We note that our analysis applies for frequency regions where only a single 
band solution exists and cu(k) is monotonic. In cases where either of these 
do not hold, multiple beams of different "rightness" may be present. 



X. HIGH INDEX VS. LOW INDEX MODULATION 

In the previous sections we focused our analysis in describing the propaga- 
tion properties of EM waves for 2D crystals with high-index contrast between 
the constituents dielectrics. We have seen anomalous refractive effects that 
include birefringence. We provided a consistent recipe based on the wave vec- 
tor diagram and band structure properties of the system, which determines 
all properties of each propagating beam such as refracted angle, phase, en- 
ergy velocity and "rightness." Since the analysis in the preceding section 
focuses on PC's with high index modulation, it is important to investigate 
the limits of validity of our theoretical analysis. We will now examine pho- 
tonic crystals with low index modulation in the context of all the properties 
that characterize the propagating beam discussed in the previous sections. 

We consider a 2D photonic crystal lattice that consists of dielectric rods 
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in air with radius 0.35 a in triangular arrangement for the H-polarization 
case. We let the dielectric constant of the rods vary, starting from the value 
of 1.05 (value close to the dielectric constant of air), and investigate the 
photonic crystal's response as the dielectric constant of the rods increases. 
Unavoidably for the low index contrast cases it is not possible to isolate 
cases where there is only one band solution for a certain frequency. For the 
five different cases that we examine in our comparative analysis, we choose 
the operation frequency to lie approximately in the middle of the spectrum 
that corresponds to the 2 nd and 3 rd band. The cases that we analyze are 
the following, a) with dielectric constant e = 1.05, b) with e = 1.2, c) with 
e = 1.5, d) with e = 2.0, and e) with e = 5.0. The operation operation 
f = fa/c is 0.80, 0.78, 0.75, 0.70, and 0.54 for the cases (a) through (e), 
respectively. We first consider the case with TK as the symmetry direction 
of the interface. The fields from the FDTD simulation for oblique incidence 
with angle 8° with the surface normal are shown in Fig. 13. 

We notice that in case (a) the wave enters essentially undisturbed inside 
the PC with the angle of propagation pretty much the same as the angle 
of incidence. As the index contrast increases, the propagating angle still 
remains close to 8°, but "wiggly" features start to appear in the phase fronts. 
Refraction angle remains positive. For index contrast 5.0 (case (e)) we see a 
beam in the negative direction. We have plotted the EFS for all five cases in 
the first zone in Fig. 14. In Fig. 15 we show the band structure for the two 
limiting cases ((a) and (e)). The solid lines correspond to the second band, 
while the dashed lines correspond to the third band. In Fig. 14 the EFS that 
are closed curves belong to the second band, while the ones that are broken 
with a six-fold symmetry belong to the third band. Notice that the second 
band EFS are highly anisotropic for ) and become more and more 

isotropic as the index contrast increases. They become eventually "almost" 
circular for case (e). If we check the FB wave components in expression (4), 
we see that for cases (a) through (c) there is mainly one predominant FB 
wave coefficient. Mixing starts to appear in case (d) and becomes stronger 
in case (e). 

Suppose that we could describe our periodic system as an effective medium 
that has a dielectric constant given by the average value of its components 
and, therefore, an effective index n e //. 

n eff = y/ e *f+(l-f) (40) 
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with f being the filling ration and e the rod dielectric constant. We note at 
this point that an ra e // corresponding to Maxwell- Garnett theory [3H] yields 
a similar propagating angle to the one predicted by the n e // of Eq. (40). 
The field solution will then be a plane wave, therefore: 

H(r,t) ~ Ae^-^. (41) 

with K having the value that corresponds to a homogeneous medium with 
refractive index n e ff i.e., K = n e ffu/c. 

Alternatively, one can consider the system as a periodic medium, with 
the field solution an FB wave given by expression (4). For cases (a) through 
(c) in the sum there is only one predominant term. Therefore, 

H(r, t) = Ae i(k P~ d - r - wt) + 0(e 2 ). (42) 

We found that k pred =k + G , k in the first BZ and G =47t/a/3 a -1 y, 
with y the propagation direction (TM in this case) (x represents the lateral 
direction). In Table I, we show for all cases (a) through (e), the wave vector 
K, assuming an effective medium with index given by Eq. (40) and the 
wave vector of the predominant component k pre d of the FB sum, that we 
determined from the PWE method. For cases (d) and (e) there are more 
than one plane wave components in the FB sum with significant magnitude. 
In these cases, for k pred , we chose the one with the larger amplitude. In 
Table I, we also show the angle of propagation 6, if we treat the system 
as a homogeneous system with index n e ff. Moreover, we show the angle 
of propagation 9 pre d assuming that in the PC the predominant plane wave 
component propagates by itself. We notice a good agreement between K and 
kpred- A small discrepancy, becomes larger as the index contrast increases 
and mixing becomes stronger. We observe that in the cases with almost no 
mixing ((a) through (c)), the angles and 9 pre d agree with the refracted 
angles observed in the simulations in Figs. 13(a) through (c). Nevertheless, 
in case (e) where mixing is quite strong, although and Q vre d are close, both 
are quite different from the actual refracted angle seen in the simulation (Fig. 
13(e)). The refracted angle in case (e) is negative. This can be understood if 
we look at the energy velocity expressions in Sec VII. For cases (a) through 
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(c), there is only one surviving term in expression (31) for Gi and G 2 = G . 
Therefore, the energy velocity becomes, 

c 2 

v e = (k + G ) fc.Go ^G ( k . w n,k) oc k pred . (43) 

TABLE I 





K x 


Ky 


|K| 


e 


kpredx 


kpredy 




Qpred 




(tt/o) 


(it /a) 


(tt/o) 


(°) 


(ir/a) 


(ir/a) 


(tt/o) 


(°) 


(a) 


0.223 


1.602 


1.618 


7.91 


0.223 


1.602 


1.617 


7.91 


(b) 


0.217 


1.613 


1.628 


7.66 


0.217 


1.610 


1.624 


7.68 


(c) 


0.209 


1.645 


1.658 


7.23 


0.209 


1.632 


1.645 


7.29 


(d) 


0.195 


1.671 


1.683 


6.65 


0.195 


1.638 


1.650 


6.78 


(e) 


0.150 


1.794 


1.800 


4.79 


0.150 


1.835 


1.841 


4.68 



So, when the index contrast is very low (cases (a) through (c)), only one 
component has a significant magnitude in the FB sum. Then, the direction 
of the energy velocity is very close to the direction of the predominant wave 
vector. However, as the mixing becomes stronger, the directions of the pre- 
dominant wave vector and energy velocity can be very different (case (e)). 
Notice, in the latter cases, other wave vectors contributing to the FB sum 
can have an amplitude quite close to the predominant one. Alternatively, if 
we would like to describe the system graphically, for cases (a) through (c), 
we could accomplish this in two different ways. We could treat the system 
as a homogeneous system and draw its EFS as a circle in wave vector space 
with radius equal to K — n e ffuj/c. But, we can consider also the medium 
as a periodic medium and follow the recipe of Sec. III. Both treatments in 
these cases give almost the same angle of propagation. This is in excellent 
agreement with what was observed in the FDTD simulation. 

One might be tempted to describe a photonic crystal medium for cases 
with a low index contrast as a homogeneous medium with an index given by 
Eq. (40). However, the results we present in the following suggest that such a 
treatment would be erroneous. In fact, consider the five different cases of Fig. 
13. In this case we take the same angle of incidence, but choose the interface 
along TM and, therefore, the propagation direction y along TK. We present 
the results in Fig. 16. Contrary to one's expectations for a homogeneous 
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medium with index n e //, even for an index contrast as low as 1.2:1.0 (case 
b), we see three distinct refracted beams. These beams have propagating 
angles in excellent agreement with the predictions of a wave vector type of 
analysis in the repeated zone scheme. Therefore we can infer that the wave 
"sees" the periodicity of the medium even when the index contrast is low. 
The periodicity may introduce multiple bands for a certain frequency (see 
Fig. 15). In addition, because the system is periodic, the modes repeat 
themselves in wave vector space. Our results indicate, that the periodicity 
of the system should be taken into account even for the low index contrast 
cases. An effective medium approach may give inadequate results (predicts 
a beam close to the position of the middle beam in cases of Fig. 16(b) 
through (e)) or, in many cases, totally inaccurate for the high index contrast 
cases. For example, for case of Fig. 13(e), it predicts a positive instead 
of a negative refracted beam. For infinitesimally small index modulation, 
however, a wave vector analysis may predict additional beams. In such a 
case the periodicity and band folding are a mere artificiality and the medium 
is essentially homogeneous and should be treated as such. 

We notice also, that in Fig. 16, for cases (b) through (c), although three 
beams are present, each has almost clear wavefronts. Actually we checked 
the field solution with the PWE method and we found that the FB wave 
describing each of these beams consists of only one predominant coefficient. 
Therefore, in these cases for each beam, we have information of how fast 
the phase given by k • r travels from point to point in space. So, in such 
cases, Yariv's picture (definition (5)) [23] is appropriate for the phase velocity. 
Therefore, it is natural to ask, when does Yariv's picture begin to fall apart? 
To answer this, we consider a similar numerical experiment as in Sec. V. 
For detector /sampling points, we use adjacent points in the numerical grid 
space. We take normal incidence and TM as the propagation direction y. 
Let yi represent the location of a point in the numerical grid space where we 
sample the field. We calculate the Fourier transformed field ratio at adjacent 
points, H (uj , yi + \) / H \u , yi) . From this ratio, we extract k which will be 



the propagation direction y. Since Ay is small in this case, we take the 
principal value of the complex logarithm in Eq. (44). For a homogeneous 




(44) 



where Ay = y i+ i — y { = \/3/62a, is the size of the numerical grid along 
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system, the extracted wave vector value would be independent of the location 
of the pair of detector points. We plot the extracted wave vector value for 
different sample points along the propagation directions. In this way, we 
can check how much the system diverges from the homogeneous case with 
increasing index contrast. For different pair of sampling points located at yi, 
we plot the extracted real k^ and imaginary part of kj of the wave vector 
k, for all five cases in Figs. 17(a) and 17(b) respectively. We find that for 
case (a) this value is almost constant. The small imaginary part present, 
is due to numerical errors. Nonetheless, as the index contrast increases, 
deviations appear from a constant value that become larger and larger. 
In addition, an increasingly large imaginary part appears, contradictory to 
the fact that the photonic crystal is an inherently lossless system. Obviously, 
a phase velocity defined in the 2D crystal according to Yariv's picture |23] 
quickly breaks down as the index contrast increases. Thus, for the large index 
contrast cases, the definition in Sec. IV for the phase velocity according to 
Notomi's picture pQ becomes appropriate. 

We now discuss the definition of "rightness" for the low index cases. This 
property was found to have a sign that coincides with the curvature of the 
band for the cases where we have a strong mixing (high index contrast). As 
we have already mentioned, the band folding is just an artificiality when the 
medium is homogeneous or when the index contrast is very low. Negative 
curvature in such a case can by no means imply the existence of a left-handed 
(backwards) beam. The question arises when is it appropriate to associate the 
sign of band curvature with the sign of the "rightness" ? All beams in Figs. 13 
(a) through (c) and Figs. 16 (a) through (c) can very well be approximated 
by a plane wave. Therefore, clearly all the observed beams in such cases are 
right-handed beams. The cases of Figs. 13(d) and 16(d) have some small 
mixing, while the cases of Fig. 13(e) and 16 (e) have a stronger mixing. 
As mentioned before, only in cases with a strong mixing does the phase 
velocity of the entire profile of the magnetic field within the Wigner-Seitz 
cell have physical meaning. Apparently, only cases with a strong mixing, 
can be candidates for left-handed behavior. In fact, we should attempt to 
discuss "rightness" only for cases where the real part k^ of the wave vector 
value (extracted using the methodology of the preceding paragraph) shows 
such large variations, that it ranges from positive to negative values. This is 
the case only for case (e) (index contrast 5.0:1.0). 

Conclusively, in the low index contrast cases, the wave vector diagram 
analysis is still valid. Even for a low index contrast in the higher bands, an 
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effective medium cannot always describe the refractive properties of the pho- 
tonic crystal. In an anisotropic homogeneous medium, at most, two beams 
would be expected but never three, as we observed in Figs. 16 (b) through 
(e). However, it is inappropriate to assign a "rightness", according to the 
"sign" of the curvature of the band structures for PC's with low index modu- 
lation. Modulation must be high enough so that no predominant plane wave 
component exists. 



XI. COMPARISON WITH THE ID LAYERED MEDIUM 

As we have mentioned before the refractive properties of the ID layered 
medium have been extensively studied (TU1 HH C2 EH El- However, there 
are significant differences between the properties of the ID layered medium 
and the two-dimensional photonic crystal we studied in this paper. In the 
two-dimensional photonic crystal, when the plane of incidence is chosen to 
be the periodic plane, the entire wave vector is confined in the first BZ. In 
contrast, in the ID periodic medium, only the component of the wave vector 
along the direction of periodicity is Bloch confined, i.e., restricted within the 
first BZ, in this case. This has several implications. 

First, the modes in k-space repeat themselves periodically in one direction 
only. If the interface is chosen along the ID periodic direction, modes from 
the higher order zones, when |k||a/7r| < 1, can never be accessed. One can 
access higher order modes for small |k||| values only when a slanted interface 
is employed. This is the method used in Refs. [TU] and [TT] to access modes 
lying outside the first BZ with small |k|||. Note, in our 2D system we can 
access higher order modes even when the interface is cut along a symmetry 
direction and small |k|||. In fact, these are propagating waves indicated with 
dotted line in the cases of Fig. 1(b) and 1(c) (type I higher order waves). 

Second, the most important implication is regarding the "rightness." In 
the 2D system, a band with negative curvature corresponds to a left-handed 
(backwards) beam. However, this is not true for a one- dimensional system. 
We have chosen x to represent the direction of the periodicity. The curva- 
ture of a certain band will then be given by dui{k X) k y )/dk x or equivalently 
v gx k x , where v gx is the component of the group velocity along the direction 
of periodicity. Following the methods used in App. II for the group velocity 
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and Sec. VII for the Poynting vector for the ID layered system, we obtain, 



< S > k = - — (k 2 A v + -^v gx k x ), (45) 

where <> refers to the spatial average within the unit cell of the time av- 
eraged quantity, u = u{k X) k y ) for the band under consideration, and if the 
fields are H-polarized 

A v = Vg,g>H g H g , = / 7](x)v 2 dx. (46) 

G,G' J 

r)(x) = and v = J2g ^G(k x , k y ) e lGx for the band under consideration. 
Since the integrand in expression (46) is positive, A n is a positive definite 
quantity [35] . 

Now, for a band with a positive curvature, v gx k x > and so < S > k > 0. 
For a band with a negative curvature v gx k x < 0. Then < S > -k < 0, only if 

tiyA-rj < -^\v gx k x \, (47) 

where k x is within the limits of the ID BZ. 

In other words, a propagating wave that corresponds to a band with 
positive curvature is always a forward wave. In contrast, a propagating 
wave that corresponds to a band with a negative curvature is backwards 
(left-handed) only when condition (47) is satisfied. Note that condition (47) 
holds, regardless of the choice of the interface, provided that x represents 
the stacking (periodic) direction and y the direction perpendicular to this. 
If we choose the interface along y and consider normal incidence, then k y = 
0. Thus, in this particular band with negative curvature yields a 

backwards wave. Furthermore, we examined this condition for with 
high index modulation, H-polarization, frequency falling in the second band, 
and interface along the x-direction. With H-polarization in this case, we 
mean that the magnetic field is perpendicular to the plane of incidence. We 
employed the PWE method [SH EE1 EE] and found that the possibility of 
left-handed behavior, is restricted only to a small fraction of frequencies of 
the second band. In addition, for the applicable frequencies, one obtains 
backwards waves only for a part of the wave vector space. So, for the ID 
layered medium, negative curvature does not necessarily imply a backwards 
beam. Each individual case should be examined with condition (47), to 
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determine the "rightness" of the propagating beam. We note at this point 
that backwards coupling jUj observed between two waveguides linked with 
ID layered medium does not necessarily imply a backwards wave. 



XII. CONCLUSIONS 

We systematically studied EM wave propagation in two-dimensional pho- 
tonic crystal structures. We based our analysis on the wave vector diagram 
formalism. We observed different cases where negative refracted beam with 
distinctly different origins are present. We confirmed the condition for single 
beam propagation does not coincide with the condition for having a single 
reflected beam in the incoming medium. For simple cases, we determined the 
conditions for single beam propagation and applicability of Snell's formula. 
We revisited the controversial topic of phase velocity and showed that in a 
photonic crystal with strong scattering present, only the wave vector inside 
the first BZ zone has physical meaning. We used the symmetrical properties 
of the photonic crystal to appropriately design a wedge experiment that can 
determine the "rightness" of a general 2D PC system (hexagonal or square). 
We studied the behavior of the PC system as the index contrast transitions 
from high to low values. We used, whenever possible, the symmetry of the 
system to determine the reflected beams by inferring a simple formula. For 
the cases in Sec. Ill the refracted beams were determined by the use of one 
primary construction line. With the rapid development of photonic crystals 
more complicated structures are now fabricated, for example 12-fold symmet- 
rical quasi-crystals |18]. In more complicated structures, or when interfaces 
are not cut along crystal symmetry directions, the wave vector diagram anal- 
ysis should be performed in its general form. All construction lines should be 
kept and all intersections should be accounted for to determine all possible 
reflected beams, as well as all possible refracted beams. Our study will help 
in the understanding of EM propagation in more complicated and/or three 
dimensional structures. In 3D structures, interesting phenomena may arise 
because of the possibility of polarization coupling. The present study will 
also help in the understanding and making of optical devices such as light 
deflection devices jlH], waveguide division multiplexers [SU] etc. 
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APPENDIX I 

a) Higher order Bragg reflected beams 

The kn m component of an order m Bragg reflected wave in the air medium 
is given by Eq. (1). In order not to have any Bragg waves for any angle of 
incidence the condition, 

.uj . 2mn . . 

\ — smd inc +- \>w/c, (48) 

must be observed Vm^0,V 9i nc e [0, n]. 
But, if 

-sm6 inc >- , (49) 



C Dstr,cut 

then 



\ U ■ n UJ . 2-k UJ . UJ 

l-sm^nc - 1 = —smB inc - < -smO inc < -, (50) 

C Vstr,cut C "str,cut C C 

i.e., a Bragg wave of order m=-l couples. Therefore, the first condition we 
must impose is 



u ■ n 27T 

-sin0 inc < , (51) 

C "str,cut 
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V \&inc\ < 71 ■ Equivalent ly, 



U 27T 

7 < ?: ( 52 ) 

C "str,cut 



must be satisfied. Assuming this condition is valid, we proceed. 

uj 2mn uj Imix uj 2n 
\-smB inc + 1 > min(\ — sm6 inc + 1) > | 1. (53) 

C "str,cut C "str,cut C "str,cut 

So, it suffices to find the frequencies to satisfy, 

> - (54) 



or equivalently, 



f<U , (55) 



with a being the lattice constant. 

No higher order reflected beams appear for any angle of incidence for 
frequencies satisfying (55). If condition (55) (or (54)) is valid, condition (52) 
is automatically satisfied. In addition, for a certain frequency satisfying (52), 
for a certain incident angle 9 inc obeying the inequality, 

- — sin 6 inc + 1-njh str cut > —, (56) 
c c 

no higher order Bragg reflected beams appear. 

b) Equivalent points in wave vector space 

The eigenvalue equation that results from expression (4) and Maxwell's 

equation is for the H polarization case, 

J> G)G ,(k + G) • (k + G')tf G ,(k) = ^H G (k). (57) 



G' 



Suppose we consider the FB wave for K = k + G Q with G Q a reciprocal 
lattice vector. Then, the eigenvalue equation becomes 

J2 VgM* + G + Go) • (k + G' + G )tf G > (K) = ^ifo(K). (58) 
g' C 
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From Eq. (58) after setting Gi = G + G Q and G 2 = G' + G Q we get 

J> Gl , G2 (k + Gi) • (k + G 2 )iJ G2 _ Go (K) = ^tf Gl _ Go (K). (59) 
g 2 c 

By comparison with the original eigenvalue equation, it is evident that 

# G _ Go (K) = # G (k). (60) 

Therefore, the time-independent part of the FB wave (Eq. (4)) for K 
becomes 

H fb ,k = e iK r #g(K) e iG r = e lk r ^ H G (K) e *( G + G °> r 



= e lkr X> G '_ Go (K) e tG ' r = e ik ' r E^c'(k) e lG '- r = H FB , k . (61) 

G' G' 

So, the Floquet-Bloch wave expressions corresponding to wave vectors k and 
K, separated by a reciprocal lattice vector G Q , are equivalent. In other 
words, k and K are equivalent points in wave vector space. 

c) Yariv's definition for phase velocity 

In the ID system the wave vector in the plane of incidence is not confined 
in the first BZ, but only the component along the periodicity. Therefore, the 
FB wave has the form (when magnetic field perpendicular to the plane of 
incidence) 

H(r,t) = e ikx e if3y ^ H G (k, u) e lGx e~ wt z. (62) 

G 

x is chosen to represent the direction of periodicity. The phase velocity 
defined in Ref. |2SJ for the FB wave, given by (62) is 

v p = - (63) 

In this expression (see Ref. |25j). k is not within the first BZ zone, but chosen 
so that |# | >\H G \ VG ^ 0. 
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d) The average of the FB wave in the Wigner-Seitz cell 



< if (r = R) > 




G-(R+r')^2 r / 




e ikr '^ff G e iGr '(iV = e~ iujt e 



< H (r = 0) >(64) 



where r' ranges within the Wigner-Seitz cell around r=R, R a Bravais lattice 
vector, and we used e lG ' R =l. 



Appendix II: Calculation of the group velocity 

For the group velocity calculation we employ the kp perturbation method 
that was first introduced for the photonic crystal by Johnson et al. [4*2] . 
Busch et al. jlSl HI] implemented this method to determine the group veloc- 
ity and "photon mass" of a two-dimensional E-polarized photonic crystal. In 
the following, we provide the basic steps of such calculation and derive simple 
expressions in terms of the FB wave's coefficients. The final expressions are 
compared with the corresponding energy velocity expressions in Sec. VII. 
a) H-polarization 

From Maxwell's equations we obtain: 



where H n k is the FB wave that corresponds to a band with index n, wave 
vector within the first BZ k and frequency uj^. The FB wave is given by 
expression (17) with t> ni k given by expression (19). Since the wave is H- 
polarized, the eigenvectors are parallel to the cylinder axis z, and k lies in 



V x (—-V x H n , k ) 
e r 




c 



(65) 
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the plane of periodicity xy. Using expression (65), (17) and (19) we get 



'e(r) ' " v ~ c 2 



V x (4-V x (e* k - r v„, k )) = ^ e ikr v n , k . (66) 



We set (rj(r) = hM, and take into account that k ■ v„ ik = and that 
(x, y) (since xy is the plane of periodicity). After manipulations, 

we get 



O k v„ jk = A„ ik v„ ik , (67) 

with 

6 k = V?7(r) x V x +iVr](r) x k x +rj(r)[ik x V x +k 2 - V 2 - i(k • V)](68) 
and 

A„, k = % (69) 

Equation (67) represents the eigenvalue equation that yields the FB wave 
states defined in Eq. (4) for a certain wave vector k and frequency uj n ^. 
Let us now consider the FB wave for the wave vector k = k + q where 
|q| << n /a . The eigenvalue equation for wave vector value k + q will 
be 

n,k+q A njk _|_qV njk _|_q. (70) 

Setting k — > k + q in Eq. (68) we get: 

6 k+q = 6 k + q-ft + 0(e 2 ) (71) 

with 

Q = -iVrj(r) - 2^(r)V + 2r?(r)k. (72) 
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So, the eigenvalue operator for k + q can be considered as the corresponding 
operator for k with a perturbation q • f2. Thus, 

. x < v n>k |ri|v nik > 

-"Vk+q — A ^,k + q • i { < <jj 

< v njk |v njk > 

Consequently, 

„ c 2 < v n>k | ri |v„ )k > 

w„,k+ q - w n>k = q -V k w = - q : . (74) 

^n,k < V njk |V njk > 

Apparently, 



c 2 < v n , k | ft |v„ jk > 
v g = ; ; ; . (75) 

2<^n,k < V nik |v njk > 

Note that < v n>k | O |v n k >= f wsc v* k O v n k cPr. The integral is over the 
two-dimensional Wigner-seitz (WS) cell. Each term of the operator f2 is 
evaluated separately in the numerator of the expression (75). In fact, using 
expression (19) for v n k we obtain 

< v n , k | - iV77(r)|v n , k >= J2 G 'v G 'H G+G >H G (76) 

G,G' 



< v„, k | - 2^(r)V|v n , k >= 2 J2 G Vg'H g+g >H g (77) 

g,g' 

< v nik |2r7(r)k|v nik >= 2 ^ k r] G > H G+G > H G (78) 

g,g' 



and 



< v nik | Vriik >= H 2 G = I. (79) 



G 



(eigenvectors are normalized to unity). Note to derive the above relations we 
used 



— *— / e lGr d 2 r = 5(G). (80) 
Awsc Jwsc 
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Finally, lumping all terms together we calculate 



c 2 



*s - 9 (2k + 2G + G') Vg ,H g H g+g , (81) 

After index manipulation (i.e., setting Gi = G + G and G 2 = G), we obtain 
v g = — V (k + G 1 )r ]Gl ^ G2 H Gl H G2 . (82) 

To obtain the above expression, we used 

^ Gl ?7Gi-G 2 #Gi#G 2 = ^ G 2 ^Gi-G 2 -^Gi-^G 2 ) (83) 
Gi,G 2 Gi,G 2 

since the expression inside the sum is symmetrical for the pair of reciprocal 
vectors Gi, G 2 . 

b) E-polarization 

From Maxwell's equations, 



V x V x E = ^e(r)E. (84) 
c 2 

The electric field satisfies Bloch's theorem, therefore, 

E = e* kr u n , k , (85) 

with 

u„, k = $> G e* G r . (86) 

V A wsc q 

Note in this case the electric field E is parallel to the cylinder axis z. This 
means ku nk = and E G = E G z. Expression (84) with (85) and (86) yields 

O k u nk = A nik u nik , (87) 
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with ma mi 

6 k = -V 2 - 2z(k • V) + k 2 , (88) 
and A nj k given by (69). After following a similar analysis as in a), we obtain 



c 2 < u„, ;U |fi|u nik > 



with 

ft = -2iV + 2k. (90) 
Expression (89), after acting f2 on u n k yields 



c 2 



£(k + G)£&. (91) 



G 
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Figure 1: Oblique incidence of EM waves at photonic crystal slabs. The PC 
system consists of dielectric rods in a hexagonal arrangement. All cases are 
with magnetic field along the cylinders (H-polarization). The parameters for 
each case (dielectric constant of rods e, radius of rods r, and dimensionless 
frequency /) are: a) e=12.96, r = 0.35a, /=0.58, b) e=20.0, r = 0.37a, 
7=0.425, c) e=12.96, r = 0.35a, 7=0.535, and finally in d) e=7.0, r = 0.35a 
, /=0.81. Note that f—fa/c = a/A, with a the lattice constant and c the 
velocity of light and A the wavelength of light in vacuum. The solid arrows 
indicate the transmitted, while the dotted black arrows indicate higher order 
beams inside the PC. 
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Figure 2: Wave vector diagram for the case of Fig. 1(a). The equifrequency 
surfaces are plotted (black solid lines) in the repeated zone scheme. The 
red solid circle represents the equifrequency surface for the air (incoming) 
medium. The green dot-dashed line is the primary (zeroth order) construc- 
tion line. The additional turquoise lines represent higher order construction 
lines. All intersections are indicated. The blue vector represents the fun- 
damental wave vector of the FB wave that corresponds to a causal signal. 
The respective energy velocity that coincides with the propagating signals 
direction is shown as the orange vector. In the insert the corresponding wave 
vector analysis in the first zone is shown. Clearly, an analysis in the first zone 
is sufficient in this case. For this cut G ox ,y = ^/dx,y with a x = a (lattice 
constant) and a y = v3a. A general reciprocal lattice vector is (2nl+n2) G ox 
+ n2 G oy , with nl, n2 integers. 
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- ^ refracted beams 

G ox 

Figure 3: Wave vector diagram for the case of Fig. 1(b). The equifrequency 
surfaces are plotted (black solid lines) in the repeated zone scheme. The 
red solid circle represents the equifrequency surface for the air (incoming) 
medium. The green dot-dashed line is the primary (zeroth order) construc- 
tion line. All intersections are indicated. The blue vectors represent the 
intersections that result in the first zone after the folding process, that cor- 
respond to causal signal (shown as the orange vectors). In the insert, the 
corresponding wave vector analysis done in the first zone is shown,. Thus, the 
latter fails to predict the second refracted beam. For this cut G ox y = 2n/a x ^ y 
with a x = y/3a (lattice constant) and a y = a. A general reciprocal lattice 
vector is nl G ox + (nl+2n2) G oy , with nl, n2 integers. 
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Figure 4: Same as Fig. 3 but for the case of Fig. 1(c). 
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Figure 5: Wave vector diagram for the case of Fig. 1(d) in the extended zone 
scheme. For this cut G OXiy = 27r/a X) y with a x = a (lattice constant) and 
a y = \^3a. A general reciprocal lattice vector is (2nl+n2) G ox + n2 G oy , 
with nl, n2 integers. Everything else is the same as in the previous figures. 
Notice that in this case a wave vector type of analysis in the first zone (see 
insert) predicts no propagating signal. 
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Figure 6: Oblique incidence at the photonic crystal slab with e = 60., r = 
0.37a, for frequency / = 0.275 that is below the Bragg condition for no 
additional reflected beams for any angle of incidence. Notice that despite 
the presence of only one reflected beam, there are two propagating beams 
indicated with the black solid and dotted arrows, respectively. 
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Figure 7: The open circles represent the imaginary part of the Fourier trans- 
formed magnetic field, sampled in time at different points along the propa- 
gation direction y^. The sampling points are separated by one period along 
the propagation direction which is the TM direction. PC has rods with 
e = 12.96 and radius 0.35 a, and the magnetic field lies along the cylinders 
(H-polarization). In both cases, the corresponding equifrequency surfaces 
are almost isotropic. Top panel (a) is for /=0.58 that belongs to a band 
with negative curvature. The bottom panel is for /=0.48 that belongs to a 
band with positive curvature. The solid lines are oc sin(%j), where k is the 
real part of the numerically calculated wave vector. Thus, k=-2.16 a' 1 for 
case (a) and 1.35 a -1 case (b). The stars represent the imaginary part of the 
Fourier transformed magnetic field for points in the neighborhood of y^ = 7b 
with b = ^/3- 
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Figure 8: Magnitude of the group velocity for cases with almost isotropic 
EFS for a photonic crystal of rods with e = 12.96 and radius 0.35, for re- 
polarization. In a) the results of the 5 th band (negative curvature band) are 
shown. In b) the results of the 4 th band (positive curvature band) are shown. 
The solid and dot-dashed lines represent the results from the k • p pertur- 
bation method for signal along the YK and TM direction, respectively. The 
solid lines with circles represent the results obtained when considering the 
system having an effective phase index n p . The index is calculated from the 
band structure (EFS surfaces) and is frequency dependent. Agreement be- 
tween the two results is excellent close to the band edge. Since the anisotropy 
increases as one moves away from the band edge, so does the discrepancy be- 
tween the two values. 
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Figure 9: Wave vector diagrams for the wedge experiment that corresponds 
to Fig. 1(b). The black solid lines represent the EFS in the PC medium. 
The red circle is the EFS in air. The blue vector represents the wave vector 
inside the PC medium. The black solid arrow represents the causal direction 
of the energy flow inside the PC. The turquoise lines indicate the directions 
of the first and wedged interfaces, both chosen along Ti^. The normal to 
the wedged interface is indicated with the dot-dashed line. The blue solid 
line indicates the construction line at the wedged interface. The green arrow 
represents the wave vector and direction of the outgoing beam (in air). 
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Figure 10: The wedge experiment for the case of Fig. 1(b). The operation 
frequency is f — 0.425. A large positive outgoing angle is seen as predicted 
by the wave vector diagram analysis. The second beam is along the normal 
to the wedge direction and is due to multi-reflections in the upper part of 
the wedge. 
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Figure 11: The wave vectors inside the PC wedge (area 2), and the outgoing 
wave vectors (area 3) for a case of anisotropic EFS. In (a) the wave vectors 
in area 2 are drawn assuming the PC is right-handed and in (b) assuming 
the PC is left-handed. Three different beams couple into the PC with three 
different wave vectors (k t , kdi, kd2)- The red circle represents the EFS in air 
for the relevant frequency (/ = 0.50). The bold and dotted blue and green 
lines perpendicular to 12 represent the different k\\ values we obtain from the 
careful study of the wave vector diagram in the repeated zone scheme. We 
have four different outgoing beams. The position of the two is very close to 
the position of the remaining two. From the figure it becomes clear that the 
location of the beam with the larger angle coincides with the sign of v e • k 
inside the PC wedge and so determines the "rightness." 
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Figure 12: Wedge simulation for the case of Fig 11. The outgoing beam 
with the larger angle is the negative hemisphere. Therefore, the system is 
left-handed in this case. 63 




Figure 13: Refraction at oblique incidence with angle of 8° at a PC lattice of 
rods with radius 0.35 a for H-polarization. We consider a dielectric constant 
equal to 1.05, 1.2, 1.5, 2.0 and 5.0 for cases (a) through (e) respectively. 
Positive refraction is seen in all cases except in case (e) where index contrast 
is high. 64 




Figure 14: Equifrequency surfaces for 5 values of index contrast (1.05, 1.2, 
1.5, 2.0 and 5.0 to 1.0 respectively). The operation frequency is chosen 
around the middle of the second and third band. The closed curves are the 
EFS that correspond to the second band, while the open concave-like curves 
are those that correspond to the third band. The EFS are drawn in the first 
quadrant only, so that we are able to see more detail. The third band EFS 
remain anisotropic with increasing index contrast. The ones that correspond 
to the second band become increasingly isotropic and finally almost circular 
for an index contrast as high as 5.0. In the right panel the EFS are drawn 
in the 27r-space for the limiting cases with e— 1.05 (dotted lines) and e=5.0. 
(dashed lines). 
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Figure 15: The band structure for the two limiting cases with e— 1.05 and 
e=5.0. The operation frequency is chosen to be around approximately the 
middle of the second and third band, and is indicated in the figure with the 
dotted line. 
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Figure 16: Refraction at oblique incidence with an angle of 8 for the same 
cases as in Fig. 13, but with YM as the symmetry direction of the interface. 
Even for an index contrast as low as 1.2, one can see three distinct beams 
(although two of them are faint in magnitude). 
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Figure 17: (a) Real and (b) imaginary part of the wave vector for the PC 
lattices of Figs. 13-16. The wave vector is calculated from the numerical 
FDTD field patterns at adjacent points yi and yi + Ay (Ay = \^3/G2 a). 
We have taken normal incidence along TM and assumed in the wave vector 
extraction that one plane wave component dominates the propagation. 
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